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Abstract 

The elliptic Monge-Ampere equation is a fully nonlinear Partial Differential Equation 
which originated in geometric surface theory, and has been applied in dynamic meteorol- 
ogy, elasticity, geometric optics, image processing and image registration. Solutions can 
be singular, in which case standard numerical approaches fail. 

In this article we build a finite difference solver for the Monge-Ampere equation, 
which converges even for singular solutions. Regularity results are used to select a priori 
between a stable, provably convergent monotone discretization and an accurate finite 
difference discretization in different regions of the computational domain. This allows 
singular solutions to be computed using a stable method, and regular solutions to be com- 
puted more accurately. The resulting nonlinear equations are then solved by Newton's 
method. 

Computational results in two and three dimensions validate the claims of accuracy and 
solution speed. A computational example is presented which demonstrates the necessity 
of the use of the monotone scheme near singularities. 

Keywords: Fully Nonlinear Elliptic Partial Differential Equations, Monge Ampere 
equations, Nonlinear Finite Difference Methods, Viscosity Solutions, Monotone 
Schemes, Convexity Constraints 



1. Introduction 

In this article we build a finite difference solver for the Monge-Ampere equation, which 
converges even for singular solutions. Regularity results are used to select a priori be- 
tween two discretizations in different regions of the computational domain. Near possible 
singularities, a stable, provably convergent monotone discretization is used. Elsewhere a 
more accurate discretization is used. This allows singular solutions to be computed using 
a stable method, and regular solutions to be computed more accurately. The resulting 
nonlinear equations are then solved by Newton's method, which is fast, ©(M^-^), where 
M is the number of data points, independent of the regularity of the solution. 
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1.1. The setting for equation 

The Monge- Ampere equation is a fully nonlinear Partial Differential Equation (PDE). 

<let{D'^u{x)) = /(x), for x in Q.. (MA) 

The Monge- Ampere operator, det(Z)^M), is the determinant of the Hessian of the function 
u. The equation is augmented by the convexity constraint 

u is convex, (C) 

which is necessary for the equation to be elliptic. The convexity constraint is made 
explicit for emphasis: it is necessary for uniqueness of solutions and it is essential for 
numerical stability. 

While other boundary conditions appear naturally in applications, we consider the 
simplest boundary conditions: the Dirichlet problem in a convex bounded subset C M'' 
with boundary dVt, 

u{x) = g{x), for x on dH-. (D) 



nder suitable assumptions on the domain and the functions f(x),g{x), recalled in sub- 
there exist unique classical (C^) solutions to ( |MA ), ([C|). However, when 



section 2.1 



these conditions fail, solutions can be singular. For singular solutions, the correct notion 

24 of weak solutions must be used. In this case, novel discretizations and solutions methods 

25 must be used to approximate the solution. 

26 1.2. Applications 



The PDE (MA I is a geometric equation, which goes back to Monge and Ampere 
(see [T]). The equation naturally arises in geometric problems of existence and uniqueness 
of surfaces with proscribed metrics or curvatures O Early applications identified 
in [4 include dynamic meteorology, elasticity, and geometric optics El [71 [8] . For an 
application of Monge- Ampere equations in mathematical finance, see [9]. 

The Monge- Ampere equation arises as the optimality conditions for the problem of 
optimal mass transport with quadratic cost pQ (TUl E] • This application of the Monge- 
Ampere equation has been used in many areas: image registration pJJ [T31 E], mesh 
generation [151 Ell HZ] , reflector design |18j , and astrophysics (estimating the shape of 
the early universe) [T9] . 

The problem here is to find a mapping g{x) that moves the measure /ii(a;) to ^2{y) 
and minimizes the transportation cost functional 

j ^ \x - g(a;)|^ dni. 

The optimal mapping is given by g = Vu, where u satisfies the Monge- Ampere equation 

det{D'^u{x)) = ^i(a;)//i2(Vu(a;)). 

In this situation, the Dirichlet boundary condition ( |D| ) is replaced by the implicit bound- 
ary condition 

g(-) : ^ (1) 
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where the sets and ^2 are the support of the measures fJ.i,fJ.2- These boundary 
conditions are difficuh to implement numerically; we are not aware of an implementation 
using PDE methods. For many applications, both domains are squares, and a simplifying 
assumption that edges are mapped to edges allows Neumann boundary conditions to be 
used. In other applications, periodic boundary conditions are used. 

In other problems, the Monge- Ampere operator appears in an inequality constraint in 
a variational problem for optimal mappings where the cost is no longer the transportation 
cost. Here the operator has the effect of restricting the local area change on the set of 
admissible mappings, see [20] or [5T]. 

1.3. Related numerical works 

Despite the number of applications, until recently there have been few numerical 
publications devoted to solving the Monge-Ampere equation. We make a distinction 
between numerical approaches with optimal transportation type boundary conditions ([T]) 
and the standard Dirichlet boundary conditions ([d]). In the latter number of 

numerical methods have been recently proposed for the solution of the Monge-Ampere 
equation. 

An early work is j4Jj which presents a discretization which converges to the Alek- 
sandrov solution in two dimensions. Another early work by Benamou and Brenier |22j 
used a fluid mechanical approach to compute the solution to the optimal transportation 
problem. 

For the problem with Dirichlet boundary conditions which is treated here, a series of 
papers have recently appeared by two groups of authors. Dean and Glowinski 
and Feng and Neilan, [261 EZ] • The methods introduced by these authors perform best 
in the regular case and can break down in the singular case. See [28] a more complete 
discussion. 

We also mention the works [29 , in the periodic case, and [15 for applications to 
mappings. The method of jBO] treats the problem of periodic boundary conditions in 
odd dimensional space. 

1.4. Numerical challenges 

When the conditions for regularity are satisfied, classical solutions can be approx- 
imated successfully using a range of standard techniques (see, for example works such 
as [231 1211 US], and [2S1I27]). However, for singular solutions, standard numerical methods 
break down: either by becoming unstable, poorly conditioned, or by selecting the wrong 
(non-convex) solution. 

Weak solutions 

For singular solutions, the appropriate notion of weak (viscosity or Aleksandrov) 
solutions must be used. Numerical methods have been developed which capture weak 
solutions: Oliker and Prussner, in [4 , presents a method which converges to the Aleksan- 
drov solution. One of the authors introduced a finite difference method which converges 
to the viscosity solution in |31j . Both of these methods were restricted to two dimensions. 
However, methods which are provably convergent may have lower accuracy or slower so- 
lution methods than other methods which are effective for regular solutions. In [ST^ we 
introduced a monotone discretization which is valid in arbitrary dimensions. A proof 



85 of convergence to viscosity solutions is provided, as well as a proof of convergence of 

86 Newton's method. 



Convexity 

The convexity constraint is necessary for both uniqueness and stability, 
ular, the equation (MA) fails to be elliptic if u is non-convex 
instabilities can arise if the convexity condition (IC]) is violated 



In partic- 



(see 


subsection 2.5 


. so 


as demonstrated in 


sub- 



[section 83 Any approximation of (MA) requires some selection principle to choose the 
convex solution. This selection principle can be built in to the discretization, as in |31j . 
or built in to the solution method, as in (5H1. 



Accuracy 

The convergent monotone schemes of [3T] and |32] use a wide stencil, and the accuracy 
of the scheme depends on the directional resolution, which depends on the width of the 



stencil. As we demonstrate below, for highly singular solutions, such as (17), the direc- 
tional resolution error can dominate. On the other hand, more accurate discretizations, 
such as standard finite differences, can be unstable for singular solutions. 



Fast solvers 

Previous work by the authors and a coauthor [55] investigated fast solvers for ( MA I . 
An explicit method was presented which was moderately fast, independent of the solution 
time. For regular solutions, a faster (by an order of magnitude) semi-implicit solution 
method was introduced (see subsection 6.2) but this method was slower (by an order of 
magnitude) on singular solutions. 



2. Analysis and weak solutions 



In this section we present regularity results and background analysis which inform 
J^he numeri cal approach taken in this work. In particular, the regularity results of [sub-[ 
[section 2T] are used to determine the discretization used in [section 5] 
P The definition of viscosity solutions and Aleksandrov solutions presented in [subsec^ 
[tion 2.2|2.3 are used to make sense of the weak solutions (15) and (17), respectively. 



2.1. Regularity 

Under the following conditions, the Monge-Ampere equation is guaranteed to have 
a unique C^'" solution Regularity results for the Monge-Ampere equation have been 
established in [531 [Ml [3S] . We refer to the book (35] for the following result. 



The domain $7 is strictly convex with boundary dO. £ C^'". 

The boundary values g € C^'"(9r2). 

The function / £ C"(il) is strictly positive. 



(2) 



116 Remark 1. In the extreme case, with f{x) — for all x e fl, the equation ( MA[ ),(|C|) 

117 reduces to the computation of the convex envelope of the boundary conditions |37l I38j . 
us In this case, solutions may not even be continuous up to the boundary and can also be 
119 non-differentiable in the interior. 
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120 Remark 2. While is it usual to perform numerical solutions on a rectagle, regularity 

121 can break down in particular convex polygons I39j . 



122 2.2. Viscosity solutions 

123 We recall the definition of viscosity solutions [ID], which are defined for the Monge- 

124 Ampere equation in f3B . 

125 Definition 1. Let u £ C{il) be convex and / > be continuous. The function u is 

126 a viscosity suhsolution (supersolution) of the Monge- Ampere equation in f2 if whenever 

127 convex (j) £ C^(f2) and xq £ ^ are such that {u — (f>){x) < {>){u — (f>){xo) for all x in a 

128 neighbourhood of Xq, then we must have 

det(7^2^(xo)) > (<)/(a;o). 

129 The function u is a viscosity solution if it is both a viscosity subsolution and supersolution. 

130 Example 1 (Viscosity solution of Monge- Ampere) . We consider an example which 

131 will later be solved numerically in two and three dimensions (sections[8]|9| . Consider ( MA I 

132 with solution and / given by 

«(x) = i((|x|-l)+f, /(x) = (l-l/|x|)+. 

133 This function, although it is not a classical solution of the Monge- Ampere equa- 

134 tion, is a viscosity solution. 

135 2.3. Aleksandrov solutions 

136 Next we turn our attention to the Aleksandrov solution, which is a more general weak 

137 solution than the viscosity solutions. Here / is generally a measure 135]. We begin by 

138 recalling the definition of the normal mapping or subdifferential of a function. 

139 Definition 2. The normal mapping (subdifferential) of a function u is the set-valued 

140 function du defined by 

Ou^Xq) — {p : u{x) > u(xq) + p ■ (x ~ Xq)}, for all x £ fl. 

141 For a set C il, we define du{V) — [J du{x). 

xev 

142 Now we want to look at a measure generated by the Monge- Ampere operator. 

143 Definition 3. Given a function u £ C{Q), the Monge-Ampere measure associated with 

144 u is defined as 

^,{v) = \du{v)\ 

145 for any set V C fl. 

146 This measure naturally leads to the notion of the generalized or Aleksandrov solution 

147 of the Monge-Ampere equation. 
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148 Definition 4. Let fi he a Borel measure defined in a convex set G K . Then the 

149 convex function u is an Aleksandrov solution of the Monge- Ampere equation 

det{D'^u) = fi 

150 if the Monge- Ampere measure associated with u is equal to the given meaure fj,. 

151 Example 2 (Aleksandrov solution). As an example, we consider the cone and the 

152 the scaled Dirac measure 

U(X) = |X| , ^{V) = TT / S{x)dx. 

Jv 

153 2.4- A PDE for convexity 
The convexity constraint ([C| is necessary for uniqueness, since without it, —u is also 

a solution of (|MA|. 



For a twice continuously differentiable function u, the convexity restriction (jCj) can be 
written as Z?^u is positive definite. Since we wish to work with less regular solutions, ([C]) 

158 can be enforced by the equation 

XiiD^u) > 0, 

159 understood in the viscosity sense |371I3H]: where Xi[D'^u] is the smallest eigenvalue of the 

160 Hessian of u. 

161 The convexity constraint can be absorbed into the operator by defining 

d 

det+{M)^l[X+ (3) 

162 where M is a symmetric matrix, with eigenvalues, Ai <...,< A„ and 

x^ = max(a;, 0). 



163 Using this notation, (MA),([C| becomes 

det+{D^u{x)) ^ fix), for x in Q (MA) + 

164 Remark 3. Notice that there is a trade off in defining ([s]): the constraint ([C]) is elimi- 

165 nated but the operator becomes non-differentiable near singular matrices. 

166 2.5. Linearization and ellipticity 

167 The linearization of the determinant is given by 

V det(Af) • N = trace (MadjN) 



168 
169 
170 



Where Madj is the adjugate [H], which is the transpose of the cofactor matrix. The 
adjugate matrix is positive definite if and only if M is positive definite. When the matrix 
M is invertible, the adjugate, Madj, satisfies 

Madj = det(M)M-i (4) 

We now apply these considerations to the linearization of the Monge- Ampere opera- 
tor. When u G we can linearize this operator as 

Vm detiD'^u) ■ V = trace ((D^w)^^,^'^^) . (5) 
6 



173 Example 3. In two dimensions we obtain 

Vm det{D^u)v = UxxVyy + UyyV^x - ^U^yVxy 

174 which is homogenous of order one in D^u. In dimension d > 2, the hnearization is 

175 homogeneous order d — 1 in D^u. 

176 The Unear operator 

L[u] = traceA(x)Z?^u 

177 is elliptic if the coefficient matrix A{x) is positive definite. 

178 Lemma 1. Let u G C^. The linearization of the Mange- Ampere operator, ([s]) is elliptic 

179 if D^u is positive definite or, equivalently, if u is (strictly) convex. 

180 Remark 4. When the function u fails to be strictly convex, the linearization can be 

181 degenerate elliptic, which affects the conditioning of the linear system When the 

182 function u is nonconvex, the linear system can be ustable. 

183 The definition of a nonlinear elliptic PDE operator generalizes the definition of linear 

184 elliptic operator. It also allows for the operators to be non-differentiable. 

185 Definition 5. Let the PDE operator F{M) be a continuous function defined on sym- 

186 metric matrices. Then F{M) is elliptic if it satisfies the monotonicity condition 

F{M) < F{N) whenever M < N, 

187 where for symmetric matrices M < N means Mx < x'^ Nx for all x. 

188 Example 4. The operator det^(M) is a non-decreasing function of the eigenvalues, so 

189 it is elliptic. 

190 3. Tiie standard finite difference discretization 

191 We begin by considering the standard finite difference discretization of the Monge- 

192 Ampere equation. For brevity, we describe the discretization in two dimensions, but this 

193 is easily generalized to higher dimensions. 

194 This discretization does not enforce the convexity condition ([C]), so it can lead to 

195 instabilities. In particular, we show in [subsection 8.1] that Newton's method can become 

196 unstable if this discretization is used. 

197 The Monge- Ampere operator has a particularly simple form in two dimensions: 

, ,^9 X d'^ud'^u ( d'^u \ ^ . ^ ^9 
ox^ oy^ \oxoy ) 

198 In two dimensions, the natural discretization of the operator is given by 

MA^\u\ = {Vxxu){Vyyu) - {Vxyuf 



{MA) 



7 



where, writing h for the spatial resolution of the grid, 



200 
201 



Remark 5. There is no reason to assume that the standard discretization converges. In 
fact, the two dimensional scheme has multiple solutions. In [28] this discretization was 
used, but the the solvers were designed to select the convex solution. 



202 4. Convergent monotone discretization 

203 The method of 31 describes a discretization of the two-dimensional Monge-Ampere 

204 equation that converges to the viscosity solution. In [35] we introduced another dis- 

205 cretization, which generalized to higher dimensions, and also converged to the viscosity 

206 solution. Both methods require the use of a wide stencil scheme, which has an additional 

207 discretization parameter, the directional resolution, explained below. 

208 In addition to being monotone, which means it is provably convergent, the latter 



method discretizes the convexified version of the equation, (MA)^ which is enough to 
ensure convergence of Newton's method. The proof of this result can be found in |32j . 

In this section we present the convergent discretization, which will be used to build 
the hybrid solver. 

Wide stencils 

When we discretize the operator on a finite difference grid, we approximate the second 
derivatives by centred finite differences (spatial discretization). In addition, we consider 
a finite number of possible directions v that lie on the grid (directional discretization). 

We consider the finite difference operator for the second directional derivative in the 
direction f, which lies on the finite difference grid. These directional derivatives are 
discretized by simply using finite differences on the grid 

'D^uUi = I „ {u{xi + vh) + u{xi - vh) - 2u{xi)) . 

Depending on the direction of the vector this may involve a wide stencil. At points 
near the boundary of the domain, some values required by the wide stencil will not be 



222 available; see Figure 1 In these cases, we use interpolation at the boundary to construct 

223 a (lower accuracy) stencil for the second directional derivative; see [31] for more details. 

224 Since the discretization considers only a finite number of directions there will be 

225 an additional term in the consistency error coming from the angular resolution dO of the 

226 stencil. This angular resolution will decrease and approach zero as the stencil width is 

227 increased. 
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(a) In the interior. 



(b) Near the boundary. 



Figure 1: Wide stencils on a two dimensional grid. 



228 ^.2. Discretization o] the convexified Monge-Ampere operator 

229 In two dimensions, the largest and smallest eigenvalues of a symmetric matrix can be 

230 represented by the variational formula 

Ai \A] = min v^Av, X^IA] ~ max v^Aiy. 

231 This formula was used in |31) to build a monotone scheme for the Monge-Ampere oper- 

232 ator, which is the product of the eigenvalues of the Hessian, by replacing the min, max 

233 over all directions, by a finite number of grid directions. 

234 In higher dimensions, the formula above does not generalize naturally. Instead, in |32j . 

235 we used another characterization, which applied to positive definite matrices. 

236 Lemma 2 (Variational characterization of the determinant). Let A be a dx d symmetric 

237 positive definite matrix with eigenvalues Xj and let V be the set of all orthonormal bases 

238 ojm.'^: 

V = ...,yd)\v,^ ± if i ^ j, W^.h = 1}. 

239 Then the determinant of A is equivalent to 

d d 
TT A, = min XivjAv^. 

240 We use Lemma [2] to characterize the determinant of the Hessian of a convex 

241 function (p in terms of second directional derivatives of (j). 

■ TT ^''^ 
{ut,...,ui)ev ^\ dvj^ 



Aei(D'^(h) = min TT v'^D'^Sv, 
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The convexified Monge- Ampere operator (MA)'^ can then be represented by simply 
enforcing positivity of the eigenvalues, which leads to the following, 




det+(D^(l)) = min TT 

To discretize the operator on a finite difference grid, restrict to the set of orthogonal 
vectors, Q, available on the given stencil. Then the convexified Monge-Ampere opera- 
tor 



{AIAy is approximated by 



MA 



{i/i... Indies , 



(MA) 



M 



Theorem 3 (Convergence to Viscosity Solution). Let the P DE (MA I have a unique 
viscosity solution. Then the solutions of the scheme {MA)'^'^ , converges to the viscosity 
solution of (MA) as h,d6,6 — > 0. 

The proof of convergence follows from verifying consistency and degenerate ellipticity 
and can be found in 1321. 



5. A hybrid discretization 

In this section we propose a hybrid discretization of the Monge-Ampere equation 
which takes advantage of the best features of each of the previous discretizations. We 



want to make use of the natural discretization (MA)^ wherever possible in order to take 
advantage of its simpli city and higher accuracy. However, we wish to use the monotone 
discretization (MA)^^ in regions where the solution may be singular in order to prop- 
erly capture the behaviour of the viscosity solution. In this way we hope to achieve the 
second-order accuracy of the simple discretization in smooth regions and the monotonic- 
ity necessary to capture the behaviour of the viscosity solution in non-smooth regions. 
We propose the following hybrid scheme. 



Discretize ( MA ) using a weighted average of the two discretizations: 

MA" = w{x)MA^ + (1 - w{x))MA'^' {MA)" 

where w : ^ [0, 1] is a weight function defined a priori from the data as follows. 

We first identify fi" which is a neighborhood of the possible singular set of u on 17, 
using conditions 

= {x&n\e< fix) <l/e}U{x edn\ dfl fiat or g{x) ^ C^'"} (6) 

where e is a small parameter, which we can take equal to h, the spatial resolution. 

Then define w{x) to be a differentiable function which is zero in an /i-neighborhood 
of fi", and which goes to 1 elsewhere. 

Remark 6. The hybrid scheme will sometimes be less accurate than the standard finite 
differences when the solution is , because it will lose some accuracy at the fiat bound- 
ary. While this might seem conservative, there are examples, (see [2H]), where the flat 
boundary causes blow up in the Hessian, so the monotone scheme is needed. 
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273 6. Explicit and semi-implicit solution methods 



274 Any discretization of ( MA I leads to a system of nonlinear equations which must be 

275 solved in order to obtain the approximate solution. 



276 6.1. Explicit solution methods for monotone schemes 

277 Using a monotone discretization F[u] of the Monge-Ampere operator, the simplest 

278 way to solve the Monge-Ampere equation is by solving the parabolic version of the 

279 equation using forward Euler. That is, we perform the iteration 

^n+l +dt{F[u'']- f). 

280 Explicit iterative methods have the advantage that they are simple to implement, but the 

281 number of iterations required suffers from the well known CFL condition (which applies 

282 in a nonlinear form to monotone discretizations, as explained in [42] ) . This approach is 

283 slow because for stability it requires a small time step dt, which depends on the spatial 

284 resolution h. The time step, which can by computed explicitly, is 0{h^). This was the 

285 approach used in |31j . 



301 
302 



6.2. A semi-implicit solution method 

The next method we discuss is a semi-implicit method, which involves solving the 
Laplace equation at each iteration. In [28' we used the identiy ([S]) to build a semi-implicit 
solution method. We showed that the method is a contraction, but the strictness of the 
contraction requires strict positivity of /. In practice, this meant that the iteration was 
fast for regular solutions, but degenerated to become slower than the explicit method 
when / was zero in large parts of the domain. 

The conditioning of the linearized equation ([5]), which affects solution time, depends 
on the strict convexity of the solution, see lemma [T] The convexity, in turn depends of 
strict positivity of /, see [subsection 2.1[ This explains why solution time of the semi- 
implicit solver depends on regularity. 

Next, we describe a generalization of the semi-implicit method to higher dimensions. 



We won't be using the method to solve ( MA I . Instead, we will use one iteration to set 
up the initial value for Newton's method. 

Begin with the following identity for the Laplacian in two dimensions. 



I AU| ^/{KU^ = ^ + + 2U^^Uyy. (7) 

So if u solves the Monge-Ampere equation, then 



|Au| = ^ul, + uly + 2uly + 2/ - V 1^'"!' + 2/ 
This leads to a semi-implicit scheme for solving the Monge-Ampere equation, used in 



Au"+i = ^2/+ |D2^"|2 (8) 



11 



To generalize this to R , we can write the Laplacian in terms of the eigenvalues of the 
Hessian: Au = Taking the d-th power, and expanding, gives the sum of 

all possible products of d eigenvalues. 

d 

(A«)'^ = d![]A, + P(Ai,...,Ad), 

i=l 

303 where -P(A) is a d-homogeneous polynomial, which we won't need explicitly. 

304 The result is the semi-implicit scheme 

Am"+^ = Vdlf + P{\i[D'^u^], . . . , Xd[D'^u^]). (9) 

305 A natural initial value for the iteration is given by the solution of 

Au° = ^/dlf. (10) 



306 7. Implementation of Newton's method 

307 To solve the discretized equation 

MA"[u] = f 

308 we use use a damped Newton iteration 



319 
320 



for some < a < 1. The damping parameter a is chosen at each step to ensure that 
the residual \\MA^{u") — /|| is decreasing. (In practice we can often take a — 1, but 
damping is sometimes needed.) 

The corrector w" solves the linear system 

(V„MA^[M"])z;"==Myl^[u"]-/. (11) 



To set up the equation (111, the Jacobian of the scheme is needed. Since the hybrid 
discretization is a weighted average of the monotone and standard discretization, and 
the weight function, w{x), is determined a priori, the Jacobian of the hybrid scheme will 
simply be a weighted average of the corresponding Jacobians. 

The Jacobian of the Monge- Ampere operator, discretized using standard finite differ- 
ences, is given by 

which is a discrete version of the linearization of Monge- Ampere ([s]) 

The Jacobian for the monotone discretization is obtained by using Danskin's Theo- 
321 rem |43j and the product rule. 

VuMA^'[u] = ;^diag I n^-.>.t« I 2?.;.; 
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(MA) 



M 



where the v* are the directions active in the minimum in 

Thus the corrector is obtained by solving the weighted average of the two hneariza- 
tions 



{w{x)WuMA^[u''] + (1 - ii;(x))V„AfA*^[u"])w" 

= w{x)MA^ [m"] + (1 - w{x))MA'^' [u"] . (13) 



323 In order for the hnear equation ( 11 ) to be well-posed, we require the coefficient matrix 

324 to be positive definite. As observed in lemma [TheorenTTI this condition can fail if the 

325 iterate u" is not strictly convex. 



7.1. Initialization of Newton's method 

Newton's method requires a good initialization for the iteration. Since we need the 
resulting linear system to be well posed it is essential that the initial iterate: (i) be 
convex, (ii) respect the boundary conditions, (iii) be close to the solution. 

In order to do this, we first: use one step of the semi-implicit scheme to obtain a 
close initial value. This amounts to solving ( 10 1 along with consistent Dirichlet boundary 
conditions ([d|. Then convexity the result, using the method of 37J. Since both the steps 
can be performed on a very coarse grid, and interpolated onto the finer grid, the cost of 
the initialization is low. 



7. 2. Preconditioning 

In degenerate examples, the PDE for u" ( 13 ) may be degenerate, which can lead 
to an ill-conditioned or singular Jacobian. To get around this problem, we regularize 
the Jacobian to make sure the linear operator is strictly negative definite; this will not 
change the fixed points of Newton's method. We accomplish this by replacing the second 
directional derivatives u,yi, with 



Here e is a small parameter. In the computations of section 8 we take e = ^4-^ x 10 



342 8. Computational results in two dimensions 

343 In this section, we summarize the results of a number of two-dimensional examples 

344 solved using the hybrid scheme described in [section 5} In particular, we are interested 

345 in comparing the computation time for Newton's method with the time required by the 

346 methods proposed in We also visualize the map generated by the gradient of the 

347 solution. 

348 These computations are performed on an x iV grid on the square [0,1]^. The 

349 monotone scheme used a 17 point stencil. 

350 When needed as part of the initialization, the convex envelope is computed on a 

351 coarse grid using the discretization described in [37j . Since the solution can be computed 

352 on a coarse grid, and interpolated, the added computational time is negligible. 



13 



8.1. Failure of Newton's method for natural finite differences 

In this section, we give an example where Newton's method breaks down when stan- 
dard finite difi^erences are used. 

We chose an example which is only singular at one point, on the boundary. Never- 
theless, this mild singularity is enough for Newton's method to break down. 

Consider the solution of (MA) in [0, 1]^, given by 



u(x) = 



/(x) = 2(2-|xf 



The gradient of the solution is unbounded on |x| = 2, in particular at the point (1,1). 
The singularity arises from the fact that / is unbounded there. 

Due to the singularity, there is an instability in Newton's method if the natural finite 
difference method is used. The iteration is initialized with the exact solution. The result 
after performing two iterations of Newton's method along with the gradient map, is 
illustrated in Figure 2 The correct computed solution is presented in Figure 3(g)|3(h)' 




y X 

(a) Solution after two iterations 




(b) Gradient map after two iterations 



Figure 2: Failure of Newton's method using standard finite differences: the solution oscillates and 
becomes non-convex. 



365 8.2. Four representative examples 

366 We have tested the hybrid method on a number of examples of varying regularity; 

367 the results are summarized in [subsection 8.4|8.3| To illustrate these results, we present 

368 more detailed results for four representative examples. 

369 Write X = (a;, y), and Xq = (.5, .5) for the center of the domain. 

370 The first example solution, which is smooth and radial, is given by 

u{^) = exp (^-^^ , /(x) ^ (1 + |x|^) exp(|x|^). (14) 

371 The second example, which is C^, is given by 

^.(x) = i((|x-xo|-0.2)+)^ /W=(l-^,^)^- (15) 
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372 The third example is the one which was used in |subsection 8.1| to demonstrate that 

373 Newton's method for standard finite differences is unstable. The solution is twice differ- 

374 entiable in the interior of the domain, but has an unbounded gradient near the boundary 

375 point (1,1). The solution is given by 



376 
377 



378 
379 



«(x) = -V2 - |x^ /(x) = 2(2 - |xh"'. (16) 



This final is example solution is the cone, which was discussed in subsection 2.3 It 
is Lipschitz continuous. 

u(x) = yix - xol, /==^ = 7r(5xo (17) 

It order to approximate the solution on a grid with spatial resolution h, using viscosity 
solutions, we approximate the measure by its average over the ball of radius /i/2, which 
gives 

'4//12 for |x-xo| < h/2, 
otherwise. 



8. 3. Visualization of solutions and gradient maps 



In Figure 3 the solutions and the gradient maps for the three representative examples 
are presented. For example (17) the gradient map is too singular to illustrate. To 



visualize the maps, the image of a Cartesian mesh under the mapping 



is shown, where {T>xU,'Dyu) is the numerical gradient of the solution of the Monge- 
Ampere equation. The image of a circle is plotted for visualization purposes, the equation 
is solved on a square. For reference, the identity mapping is also displayed. 

In each case, the maps agree with the maps obtained using the gradient of the exact 
solution. 



390 8.4^. Computation time 



The computation times for the four representative examples is presented in Table 1 



The computations time are compared to those for the Gauss-Seidel and Poisson iterations 
described in ^25]. The Newton solver is faster in terms of absolute solution time in each 
case. [Table 2| presents of order of magnitude solutions times. The order of magnitude 
solution time for Newton's method is independent of the regularity of the solutions and 
faster than both of the other methods. 



397 8. 5. Accuracy 

398 Numerical errors are presented in |Table"3l We compare the accuracy of the hybrid 

399 scheme to the standard finite difference discretization, (using the results of [55]) and to 

400 the monotone scheme which was also solved using Newton's method. 

401 We discuss each example in turn. 
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Figure 3: Solutions and mappings for the |(a)|(b)| identity map, [(c')|(d)| C^ example, |(e)|(f)| Cl example, 
and I (g) I (h) I example with blow-up. 
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(Lipschitz) Example (17l 
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Table 1: Computation times for the Newton, Poisson, and Gauss-Seidel methods for four representative 
examples. 
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Method 



C^^" (14) 



Regularity of Solution 
C^'" (flSl) and (fiel) 



C^'^ (171 



Gauss-Seidel 
Poisson 
Newton 



Moderate 
Fast 
Fast 

0(Af 1-3)) 



) 



Moderate 
(~0(Mi-9)) 
Fast-Slow 
C'(Mi-4)-blow-up 
Fast 

0(Mi-3)) 



) 



Moderate 

0(M2)) 

Slow 
0(Af2)-blow-up) 

Fast 
(~ 0(Afi-3)) 



Table 2: Order of magnitude computation time for the different solvers in terms or the regularity of 
solutions. Here M = N"^ is the total number of grid points. 



The solution (14) 



The standard finite difference schemes gives 0{h ) accuracy (see [21]). In this case, 
the hybrid scheme is slightly less accurate, because the monotone scheme is used near the 
boundary. On a strictly convex domain the hybrid scheme would reduce to the standard 
discretization and achieve the same accuracy. 

The effect diminishes as the number of grid points increases so that the number of 
interior points using the higher order scheme dominates. Accuracy approaches 0{h?) as 
the number of grid points increases. This is a definite improvement over the monotone 
scheme, which has its accuracy limited by the stencil width. 



The C solution ( 15 ) 



The accuracy is 0{h), which is similar to the standard discretization and better 
than the limited accuracy permitted by the monotone discretization with a fixed stencil 



414 width. We also look at the error at each point (see Figure 4); it is evident that the 

415 singularity around the circle is the factor that most affects the accuracy. Because of this 

416 non-smoothness, there is no reason to expect our scheme to produce the 0{h'^) accuracy 

417 that is possible on solutions. 



418 The blow-up solution ( 16 ) 



In this case, the hybrid scheme accuracy is 0{h^'^). This is better than the accuracy 
of both the standard discretization, which was 0{h'^'^) (28|, and the monotone scheme, 
which is limited by the stencil width. 



422 The cone solution (17 1 



For this singular example, the hybrid scheme is identical to the monotone scheme 
(since the right-hand side is either or very large everywhere in the domain). Con- 
sequently, the angular resolution (stencil width) limits the accuracy of solutions. We 
observed that the 17 point stencil reduced the error by an order of magnitude compared 
to the 9 point stencil. This dependence on the stencil width is also evident in the surface 



plot of error (Figure 4 1, which demonstrates that error is largest along directions that 
are not captured by the stencil. Since this solution is so singular the reduced accuracy 
is to be expected. 
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Example with blow-up (161 
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-3 
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0.07 X 10" 


-3 


361 
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N 



C°'^ (Lipschitz) Example ([17 1 
Maximum Error 



Standard 
10 X 10"^ 
8 X 10"^ 
6 X 10"^ 
4 X 10"^ 
3 X 10"^ 
2 X 10"^ 



Monotone 
3 X 10"^ 
3 X 10"^ 

3 X 10"^ 

4 X 10"^ 
4 X 10"^ 
4 X 10"^ 
4 X 10"^ 
4 X 10"^ 



Hybrid 

3 X 10"^ 
3 X 10"^ 

3 X 10"^ 

4 X 10"^ 
4 X 10"^ 
4 X 10"^ 
4 X 10"^ 
4 X 10"^ 



31 

45 

63 

89 

127 

181 

255 

361 



Table 3: Accuracy for the standard, monotone, and hybrid discretizations for four representative exam- 
ples. 
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431 9. Computational results in three dimensions. 

432 In this section, we demonstrate the speed and accuracy of the hybrid Newton's method 

433 for three dimensional problems. These computations are performed on an iV x x 

434 grid on the square [0, 1]"^. The monotone scheme used a 19 point stencil. 

435 The size of the computation was restricted by the available memory, not by solution 

436 time (the computations were performed on a recent model laptop). 

437 The solution methods of ^28j were restricted to the two-dimensional Monge-Ampere 

438 equation, so we are no longer able to compare solution times to Newton's method for 

439 these examples. 

440 As before, we provide specific results for three representative examples of varying 

441 regularity. In this section we use the notation 

X = {x,y,z) 

442 and let xg = (.5, .5, .5) be the centre of the domain. 

443 The first example is the solution given by 




z.(x)=expl ^ I, /(x) = (l + |x|^)exp(-|x|^). (18) 



444 The second example is the solution given by 



/(x) = 



.(x) = l((|x-xo|-0.2)+)^ (19) 

+ |x-xo|>0.2 
otherwise. 



446 The third example is the surface of a ball, which is differentiable in the interior of the 

447 domain, but has an unbounded gradient at the boundary. 

u{k) = -^3~\^\\ /(x) = 3(3-|x|V/2. (20) 
20 



As indicated by the results in [Table 4[ the hybrid Newton's method continues to 
perform well in three dimensions. (The fact that the solver required only one iteration 
for Example ( 19 ) was simply an artifact — for larger problems sizes more iterations were 
required. 



Example ([18 1 
N Max Error Iterations CPU Time (s) 



7 
11 
15 
21 
31 



N 

11 

15 
21 
31 



0.0151 
0.0140 
0.0129 
0.0121 
0.0111 



0.04 
0.10 
0.71 
6.72 
86.63 



Example ([19) 
Max Error Iterations CPU Time (s) 



0.0034 
0.0022 
0.0016 
0.0009 
0.0005 



0.02 
0.09 
0.22 
1.03 
17.12 



Example with Blow-up ( 20 1 
N Max Error Iterations CPU Time (s) 



7 
11 
15 
21 
31 



9.6 X 10"^ 
5.2 X 10^3 
4.6 X 10^3 
4.0 X 10^3 
2.9 X 10^3 



0.03 
0.11 
0.48 
7.42 
138.74 



Table 4: Maximum error and computation time for the hybrid Newton's method on three representative 
examples. 



452 10. Conclusions 

453 The purpose of this work was to build a fast, accurate finite difference solver for the 

454 elliptic Monge- Ampere equation. 

455 A hybrid finite difference discretization was used which selects between an accurate 

456 standard finite difference discretization and a stable (provably convergent) monotone 

457 discretization. The choice of discretization was based on known regularity results which 

458 depended on the boundary data, g, the right hand side function /, and strict convexity of 

459 the domain. Wherever the requirements on the data are not met, the hybrid discretization 

460 chooses the monotone discretization. 

461 The discretized equations were solved by Newton's method, which is fast, C'(M^ '^), 

462 where M is the number of data points, independent of the regularity of the solution. The 

463 implementation of Newton's method was significantly (orders of magnitude) faster than 

21 



464 the two other methods used for comparison. The hybrid discretization was shown to be 

465 necessary for stabihty of Newton's method: an example with a mildly singular solution 

466 showed that the standard discretization leads to instabilities. 

467 The hybrid discretization was introduced to improve the accuracy of the monotone 

468 discretization on regular solutions. This expected improvement was achieved. On regular 

469 solutions the hybrid solver was (asymptotically) as accurate as the standard finite dif- 

470 ference discretization. For one moderately singular example the hybrid solver was more 

471 accurate than standard finite differences by 0{h). 

472 The discretization and solution method used was not restricted to two dimensions. 

473 This allowed for the solution of three dimensional problems on moderate sized grids, with 

474 the problem size limited by computer memory, not solution time. 

475 In summary, the solver presented used a novel discretization in general dimensions, 



476 accompanied by a fast solution method. The resulting solver is a significant improvement 

477 over existing methods for the solution of possibly singular solutions of the elliptic Monge- 

478 Ampere equation, in terms of solution time, stability, and accuracy. 
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